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Regular polynomial interpolation and approxi- 
mation of global solutions of linear partial dif- 
ferential equations 

1^ ■ Jorg Kampen 

O' 

Abstract. We consider regular polynomial interpolation algorithms on recur- 
sively defined sets of interpolation points which approximate global solutions 
of arbitrary well-posed systems of linear partial differential equations. Con- 
vergence of the 'limit' of the recursively constructed family of polynomials 
to the solution and error estimates are obtained from a priori estimates for 
some standard classes of linear partial differential equations, i.e. elliptic and 
•^r ■ hyperbolic equations. Another variation of the algorithm allows to construct 

polynomial interpolations which preserve systems of linear partial differential 
equations at the interpolation points. We show how this can be applied in 
,S^ • order to compute higher order terms of WKB-approximations of fundamental 

"ti ' solutions of a large class of linear parabolic equations. The error estimates are 

sensitive to the regularity of the solution. Our method is compatible with re- 
cent developments for solution of higher dimensional partial differential equa- 
tions, i.e. (adaptive) sparse grids, and weighted Monte-Carlo, and has obvious 
applications to mathematical finance and physics. 
> 

I ■ 

1. Introduction 

l> 

f^ ' This work shows how multivariate interpolation techniques can be combined with 

00 . analytic information of linear partial differential equations (i.e. a priori estimates 

^^ ' and/or WKB representations of solutions) in order to design efficient and accu- 

rate numerical schemes for solving (systems) of linear partial differential equations. 
l^ ' These schemes are nothing but sequences of multivariate polynomials which are 

^ ^N , constructed recursively such that they solve a given linear system of partial differ- 

ed ' ential equations on a finite discrete set of interpolation points. However, additional 

information is needed in order to ensure that the sequence of interpolation poly- 
nomials converges to a (or, if uniqueness is proved, the) global solution of a given 
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linear system of partial differential equations. As we shall see, this information can 
be provided by a priori estimates which in turn lead us to error estimates in regular 
norms dependent on the regularity of the solution. We examine the situation in 
the case of linear elliptic equations with variable coefhcients. Another possibility 
is that (more or less) explicit representations of solutions are known which lead to 
problems which are easier to solve. A prominent example is the WKB-expansion 
which was investigated in [6J . The recursive structure of WKB coefficient functions 
and the error analysis lead us to the problem of regular polynomial approximation. 
In this introductionary Section we our method on an abstract level. 

1.1. Regular polynomial interpolation 

Since we are interested in the relationship between multivariate polynomial in- 
terpolation and approximation of solutions of partial differential equations, our 
focus will be on multivariate polynomial interpolation. However, in order to make 
basic ideas more accessible we shall describe algorithms in the univariate case first 
and then generalize to the multivariate case. It is well known that polynomial 
interpolation in the multivariate case is quite different from the univariate case 
in general. However, in our approach which aims at solving linear systems of par- 
tial differential equations or aims at supplementing certain strategies of solving 
partial differential equations many features are already present in the univariate 
framework. In order to avoid misunderstandings, we dwell a little on this point. 
Classically, the problem of multivariate interpolation can be stated as follows (cf. 

my- 

Given a set of interpolation points Q = {xi, ■ ■ ■ ,xn} and an N-dimensional 
space Pq of polynomials find, for given values j/i, • • • , un, a. unique polynomial 
f (z P such that 

(1.1) f{x,)^y,, jel,--- ,N. 

In this form it turns out that there is an intricate relation between sets of in- 
terpolation points and interpolation spaces that must be satisfied in order that 
the problem can be considered to be well-posed. Either we have to make some 
restrictions concerning the set of interpolation points Q (cf. |llj ) or we consider 
Q to be fixed and consider the problem of constructing the polynomial space Pg 
(cf.[l]). This amounts to a construction of the map 

(1.2) e ^ Fe 

with additional constraints such as minimality of degree (cf. [lll fT]) or monotonic- 
ity (cf. Ij). In this paper we are interested in interpolation algorithms with the 
following features 

• there are no essential restriction on the discrete set Q of interpolation 
points except that Q C D, where D is the domain of the function to be 
interpolated. 

• the map O — > Pq is monoton (indeed our basic algorithm is an extension 
of multivariate versions of Newton's interpolation algorithm). 
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• the algorithm can be extended to vector valued interpolation functions 
g : D C M" -^ M'' and if g satisfies a system of linear partial differential 
equations, then the interpolation polynomial p solves the same system of 
linear partial differential equations on the given set Q of interpolation 
points. 

• the algorithm is numerically stable and practical with respect to the prob- 
lem that the interpolation function / and arbitrary set of partial deriva- 
tives of / are to be interpolated simultaneously. For the application of 
higher order approximation of the fundamental solution of linear parabolic 
equations we comute accurate approximations of derivatives of smooth 
functions up to order 10 in order to obtain an approxmation of order 5 of 
the WKB-expansion of the fundamental solution. 

• the algorithm can be refined in order to solve well-posed linear systems of 
partial differential equations directly. 

• the algorithm can be combined with collocation methods in an efhcient 
way; it can be partially parallelized. 

• the algorithm allows for error estimates which depend on the regularity of 
the solution such that the algorithm is compatible with methods for higher 
dimensional problems of linear systems of partial differential equations 
such as sparse grids, adaptive sparse grids, and weighted Monte-Carlo. 

First we consider the problem of polynomial approximation p of a regular (i.e 
smooth or finitely many times differentiable) function 

(1.3) / : L> C M" ^ R 

defined on discrete subset of 8 C D where for m given linear partial differential 
operators 

(1.4) L, = J2 <(^)^- 

|a|<ft 

we require that 

(1.5) Lif{xj) — Lip{xj) for 1 < i < m 

for some finite set of points Xj € & G D. As indicated above we shall allow that 
the interpolation set Q can be constructed recursively (and, hence, extended arbi- 
trarily within the domain of the interpolation function). Investigations of specific 
instances of this problem can be found in the literature on polynomial interpola- 
tion (cf. the survey paper of [TU] for the development up to the year 2001). Note 
that other algorithms of natural interpolation of C'^-functions have been proposed 
(cf.[5] for hints at the history and further references). 

The paper is organized as follows. In Section 1.2 we introduce the partial 
differential equations for which we seek global regular interpolation polynomials 
of their global solutions. All basic types of partial differential equations, i.e. ellip- 
tic equations, parabolic equations, and hyperbolic equations are considered. While 
the basic algorithm is quite similar for each type of partial differential equation. 
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we shall see, however, that the convergence of the scheme of recursively defined 
interpolation polynomials depends on very different a priori estimates for differ- 
ent type of equations. In case of second order elliptic equations classical Schauder 
boundary estimates can be used, while in the case of hyperbolic equations en- 
ergy estimates are considered. In the case of parabolic equations we refer back to 
Safanov-Krylov estimates considered in the context of the truncation error analy- 
sis of WKB-expansions. In Section 2.1 we introduce first an extension of Newton's 
polynomial algorithm which interpolates a given function and its derivatives up to 
some given order k simultaneously. Section 2.2. describes a variation of this algo- 
rithm which interpolates a given function such that a given set of partial differential 
equations is preserved. Section 3 discusses the extension to the multvariate case. In 
Section 4 we refine the algorithm and construct polynomials which satisfy a given 
linear (i.g. partial) differential equation on a given set of interpolation points, i.e. 
there is no given function to be interpolated. In Section 5 we consider refinements 
which show how polynomials constructed on disjoint sets of interpolation points 
can be synthesized in order to get one polynomial which interpolates on the union 
of sets of interpolation points. Naturally, parallelization is consideredin this con- 
text. In Section 6 we show how a priori estimates of elliptic equations (standard 
Schauder boundary estimates) and hyperbolic equations (energy estimates) lead 
to convergent schemes implied by error estimates. Section 7 discusses a special 
use of regular polynomial interpolation for parabolic equations where the global 
solution is given in the form of a WKB- expansion. Section 8 provides a numerical 
example of global regular polynomial interpolation of a locally analytic function 
up to the third derivative. In Section 9 we provide a summary and give an outlook 
on current research and research in the near future. Before we start with the de- 
scription of the algorithm, we state the typical linear partial differential equations 
and indicate the different types of approximations and error estimates which we 
aim at. 

1.2. Regular interpolation and partial differential equations 

We consider the three standard types of linear partial differential equations, namely 
elliptic equations, parabolic equations, and hyperbolic equations, and exemplify 
different types of application and extension. 

• The most popular examples of elliptic partial differential equations arc of 
the second order form, i.e. 



to be solved on a domain Q, C M" with the boundary condition 

an 



(1.7) 
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for some function / : dQ -^ M which is usually assumed to be Lipschitz 
continuous at least. Here, Ujk are (at least) measurable coefficient functions 
satisfying for some constant c, and ellipticity means that 

(1.8) z_,'^jk{x)(,i^j > c > (uniformly in x). 

jk 

We construct an extension of the polynomial interpolation algorithm which 
produces a multivariate polynomial solving this elliptic equation on an ar- 
bitrary grid of interpolation points. In order to obtain error estimates b 
standard boundary Schauder estimates in this paper we shall make some 
regularity assumptions. We derive convergence of the family of multivari- 
ate polynomials constructed by our our interpolation scheme to the global 
solution of the linear elliptic equation on a bounded domain and we derive 
error estimates from a priori estimates. 
• Parabolic equations of the form 

(L9) |-i-»- 

on D := 17 X (0,T), {n C M", with 
(1.10) u(0, x) = Sy{x) ■- 5{x -y), ye M", 

where 6 is the Dirac delta distribution, and where 

1 fill //7/ 

ij •' i 

is an elliptic operator. The solution of this equation is called fundamental 
solution, because solutions of standard parabolic initial-value boundary 
problems can be represented by convolution integrals of data functions 
with the fundamental solution. The standard assumptions for such a fun- 
damental solution to exist are 

(A) The operator L is uniformly parabolic in M", i.e. there exists < A < 
A < oo such that for aU ^ e R" \ {0} 

n 

0<A|eP< 5]a,,(x)C.e, <A|^p. 

(B) The coefficients of L are bounded functions in K." which are uniformly 
Holder continuous of exponent a (a G (0, 1)). 

If some regularity assumptions on the coefficients hold in addition, then 
it can be shown that the fundamental solution p is of the form 

(1.12) p{t, X, y) ^ -^ exp - "^ ^"^l '^^ + ^ c^ (a:, y)t' 

v27rt \ ^t ^^^ 
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with some regular coefficient functions d? and Cfe. We sliall show how our 
regular polynomial interpolation algorithm can be used to compute the 
fundamental solution in terms of this representation. 

Remark 1.1. The algorithm designed in the case of elliptic equations can 
be applied to the parabolic case directly, of course. However, it turns out 
that the convergence is better if the special representation (|1.12p is used. 

• As an example of a hyperbolic equation we consider an equation of the 
form 

(1.13) Lu = /inf7, 

where 

ij ■^ i •' 

and {hij) is a symmetric matrix of signature (n, 1), if dimf7 = n + 1. We 
assume that some O <Z Vl is bounded by two spacelike surfaces Ei and 
Eg and swept out by a family of spacelike surfaces Ee(s)- We assume the 
initial conditions 

(1.15) u — 9 and du — lo 

where (7 is a function on O and fi is a 1-form. 



2. Interpolation algorithm (univariate case) 

We start with the description of the algorithm which produces polynomials which 
satisfy some given requirements on interpolation points. Our starting point is an 
extension of Newton's polynomial interpolation method such that the interpola- 
tion polynomial and its derivatives up to a given order k (an integer) equal a 
given function and its derivatives up to order k at the interpolation points. For 
simplicity of representation and since the essential features of the algorithm can 
be demonstrated for one dimensional functions, we describe our ideas first in the 
univariate case and then generalize to the multivariate case in the next section. 

2.1. Extension of Newton's method 

Let us recall the Newtonian interpolation for an univariate function 

(2.1) / : [a,6] cR-^M. 

Given a discrete set of interpolation points D — {xq, xi- ■ ■ , xn} C [a, b] we want 
to construct a polynomial 

p : [a, 6] C M ^ M such that 
(2.2) 

f{xi) = p{xi) for aU Xi e D. 
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The idea of the basic Newton interpolation algorithm is that instead of looking for 
some polynomial of form X]i=i ^i^^ ^'^^ some constants bi we may write 



(2.3) 



N 

E 

1=0 



ai'^iix) 



with 

(2.4) <i>o{x) ^ 1 and $;(x) = U-{^o{x - x^) for I > 1. 

In order to determine oq, • • • a^ we then may solve the system 



(2.5) Roa :-- 












</)i(xi) 








(t'lix2) 


02(2:2) ••• 






1 0i(xAr) (f)2{xN) 



0Ar(x7v) 



flo 
fli 
0-2 

ON 



fi^o) 

.f{xN) 



This leads to an i^-approximation of the function / similar to the Gaussian al- 
gorithm. Note however, that the matrix Rq is a lower diagonal. Hence the linear 
system can be solved easily. Moreover the matrix condition number is much better 
than that of the Vandermonde matrix used in the classical Gaussian interpolation. 
We extend this idea to a C'^-norm interpolation, i.e. we design an algorithm that 
approximates / up to the fc-th derivative, i.e. we construct a polynomial 

g : [a, 6] C K ^ M such that 
(2.6) 

f^^Hx^) = q^^\x,) for aU x, G L> and aU I < fc, 

where for a function g : [a, 6] C M ^- M y*-'-* denotes the derivative of order I while 
g = g'^ . We consider the polynomial 



(2.7) 



(Af+l)(fc+l)-l 

E 

m— 



am^m.k{.x) 



where 

(2.8) $™,fc(x) = (x - x„ div(fe+i)3 



mod(fc+i)jjmdiv(fc+i)-i. _ . 



fc+i 



where, by convention, we understand 

(2.9) nrio(a; - xif+' := 1. 

For simplicity of notation we sometimes use the abbreviations 

(2.10) p{m) = ?Tidiv(fc + 1) and q{m) ~ mmod{k + 1). 
Next we define 



(2.11) 



.(0 



K!A-) 



d 



^m,k{x), 
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and for each fc > 1 the linear system 



(2.12) 



Rk 



flo 






ai 






0-2 


=^ 




a(k+i){N+i)-i 




. f'''"Hxf^k+l){N 



where i?fc is a (iV + l)(fc+l) x (iV+l)(fc+l)-matrix determined by (fc+1) x (fc+1) 



matrices A^™ as folfows: 



(2.13) 



Rk 



AT 

Af 



AV 
Af 



j^NQ J^Nl 



■Z'fe 
■Z'fe 

Af 



AT 



Zk 
Zk 
'Z'k 



^f^ 



'Z'k 
Zk 
Zk 



j^NN 



where Zk is the (fc + 1) x (fc + 1) matrix with entries, and 

Ak \^j ) 



(2.14) 

with 
(2.15) 



A'^ ■■- 



Note that 



(2.16) 



Al 



'S(fc+i)p(i),fc(a;j) $(fe+i)p(i)+i,fc(a;j) $rfc+i)p(i)+2,fc(a;j) 



^(fe+i)p(i),fcW) 



$ 



(fc) 



*ffc+i)p(i)+i,fc(^i) 
*(fc+i)p(i)+i,fc(^i) 



,(fc) 



^ffc+l)p(i)+2,fcW) 
*(fc+l)p(i)+2,feV^j) 



,(fc) 



*ffe+i)p(i)+fc,fc(a;j) 

^(k+l)p(i)+k,ky-^3J 
^(fc+l)p(i)+fc,fcV-^j; 



(fc) 



(fc+l)p(i),fcV-^J/ ^(fc+l)p(i) + l,fc*>'^J'' ^(fc+l)p(i)+2,fcV-^J/ ^(fc+l)p(i)+fc,fcV-^J/ 






1 







1 

2 


•• 
•• 
•• 


• 

• 

• 








•• 


• fc! 



This leads to a system which can be solved row by row. It is therefore very easy 
to implement and numerically well-conditioned. 

Remark 2.1. In order to avoid large entries in the matrices Ajj" one may consider 
basis functions of form jr^/^ i ^-i (j-i /j, but we do not deal with the peculiar niceties 
of computation here. 
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2.2. Interpolation preserving linear systems of differential equations 

The preceding algorithm can be adapted it in order to construct a polynomial 
approximation p of / where the k differential operators 

(2.17) LJix) = Y. ^]i^)^f{^),i^h---,k 

are preserved on a discrete set of points O = \xq^ • • • , x^r} in the sense that 

(2.18) Lif{xj) = Lip(xj) for Xj e 6. 

At this point the linear system of the operators {Li\l < i < fc} is quite arbitrary; 
we just assume that the operators are defined pointwise, i.e. x —^ a^{x) are classical 
functions which can be evaluated pointwise (at least on the set of interpolation 
points). Note that we do not ask about convergence of a family of interpolation 
polynomials to at this point. There are several possibilities to extend our preceding 
algorithm. One is the following. Let 

(2.19) Q, := {j\a) ^ O} 
and define 



(2-20) LT= y: <^d^.- 

We start with 

(2.21) Qi = {iii,--- ,^iri}, 
and assume that 

(2.22) ill < • •• < iirt 

We consider first the interpolation point .tq and start with the following ansatz for 
the interpolation polynomial 

(2.23) Pi^ix)^ J2 bl^i^-^^oT'^- 

We assume f{xo) — Pio{xo) = w.l.o.g. ; we shall see later how we interpolate 
values of / different from zero at the other interpolation points xi,- ■ ■ ,xj^. First 
we apply the operator 

(2.24) L\^^alix)^ 
to / and pio at xq. This leads to 

(2.25) ^^lbl^ = «i^(,„)^(,„) ^ b}^ ^ ^<^(-o) 
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Inductively we assume that the coefficients b]'^ have been defined up to the index 
ifn for some m < ri and that the operator L*j^™ has been defined accordingly. We 
apply the operator 

(2.26) L\"^+' = V{- + a] ^, {x)4^^ 

to / and piQ at xq- For an integer s with m + I < s < ri define 

s 

(2.27) p^o(^)-E^'"(^-^o)'^- 

j=i 

Then we have 

i'r^Vio(a;o) = i'rPio(a:o) +aL+i(^o)£^Pio(a;o) = 
(2.28) 

This gives b}'^ . Next inductively assume that an interpolation polynomial pik has 
been constructed which interpolates / on the set of interpolation points {xq, ■ ■ ■ ,Xk} 
for some positive integer k with k < N subject to the condition 

(2.29) Lifix,) = pikix,) for 1 < ^ < fc. 

First we extend that polynomial in order to interpolate / at the point Xk+i- We 
consider the ansatz 

(2.30) pl^^^^^{x)^p^k{x)+bl^''+''^nt^{x-xiYK 
We then get b^ ' from the equation 

(2-31) p?(fc+i)(xfe+i) == I{xk+i). 

The ansatz for Pi[k+i) (i-e. the interpolation polynomial which preserves Lif on 
the set of interpolation points {xi, • • • , Xk+i}) is 

and the determination of coefficient constants b- is similar to the procedure 

for the interpolation point xq described above. Proceeding inductively, we are lead 
to the polynomial pi which interpolates / at the interpolation points of G = 
{xq, • • • , xn} such that 

(2.33) Lipi{xj) ~ Lif{xj) for all Xj e Q. 

Finally assuming that for some integer s < k the polynomial ps satisfies the 
condition that 

Ps{xj) = f{xj) for Xj G Q 
(2.34) 

LiPs{xj) — Lif{xj) for Xj G Q and i < s, 



(2.32) Piik+i){x) ^ p"^k+i)i^) + J2 bll'+'\x~Xk+ir^nl,{x^xir+' 
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it is clear that we only need to consider the reduced operator 

(2.35) L.+i^ Y. <'(^)'^^' 

ijieQs+i\u|^iQi 

and proceed analogously. 



dx^^ 



3. Extension to the multivariate case 

Next we consider generalizations to the multivariate case. There are several pos- 
sibilities but the most simple seems to be the following. First we formulate the 
problem in a way that will turn out to be useful in the context of polynomial 
interpolation of global solutions of linear systems of partial differential equations. 
In its most simple form it is a form of multivariate Newton interpolation: given a 
function 

(3.1) / : S* C R" ^ K 

we want to construct a polynomial 

p: S CR" ^M. such that 
(3.2) 

f{xi) — p{xi) for all Xi (z D <Z S, 

where D — {xq, xi, • • • , Xn} is some discrete sets of points in M" whose coordinates 
will be denoted by superscript indices as xj, j — I,- ■ ■ ,n. This is done then by 
recursive definition of polynomials po,pi, ■ ■ ■ . First, define 

(3.3) po{x) = f{xo). 
Next, ansatz and equation 

(3.4) p,ix) = f{xo) + a,U^^,ix' - x^) = /(xi) 
leads to the determination of pi by 

(3.5) ai = l,^:'']'!^'"l\ 

Next assume that pojPij ■ ■ ■ ^Vq have been defined. Then ansatz and equation 

(3.6) Pq+\(xq^i) =p(a;g+i) +aq+in^^onf^i(a:;* -4) = /(xg+i) 
leads to the determination of p^+i by 

(1 n\ „ — /(g:i;+l)-Pq(a:q+l) 
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3.1. Extension of Newton's method 

Next wc extend a multivariate version of Newton's method, i.e. we design an 
algorithm that approximates / up to the /3-th derivative (/3 — (/3i, • • • , /3n) being 
some multiindex) where we construct a polynomial 



q:S C 



such that 



(3.8) 



^(x,) = ^{x,) for &l\x,eDCS and aU 7 < /?. 



where (3 is given (i.e. the multivariate substitute for k in the univariate case de- 
scribed above), and ordering is in the following sense: 

Definition 3.1. Let x" and x'^ be monomials in M [xi, • • • , a;„]. We say that x" > x^ 
( lexicographical order) if ^^ a* > ^^ /3* or ^- a* = ^^ /3% and in the difference 
a — /3 G Z" the left-most non zero entity is positive. 



Now, let Qfo, ai 



an enumeration of multiindices with respect to 



this ordering. We define a sequence of polynomials pa„,Pai, ■ ■ ■ ,Pa^, ■ ■ ■ recur- 
sively. First, let 



(3.9) 



7</3 



If Pao, • • • ,Pa„_i have been defined, then we define 
(3.10) 

Z.^7</3 '^"m-l7l-'-i=l(^ ~^a™_i) "^j=Q "^i=\\^ '^'^ctj) 

This leads to a linear system to be solved for a vector (aqq, • • • , aaj^fj) of length 
(iV+l)(E,/3' + l) 



(3.11) 



with 



(3.12) 



R6 



Rt 



O-aoP 



O-unP 



JK^aa) 



.f{Xai) 



f^'Hx^, 



Af 


Zp 


Zp 


Zp 


Af 


411 

^/3 


Zp 


Zp 


Af 


421 


431 

^/3 


Zp 



aNH aNI aN2 aNZ 

Ap Ap A^ /i^ 



Zp 

Zf3 
Zf3 



ANN 
^/3 
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We abbreviate J2 f^ — J2iil^^ + 1) ^^^^ defining p{m) = m -f- ^ /? we have 
(3.13) 



$ 



Y.I3p{i),l3y-^3) ^E/3pW+ft,/3'^ J'' ^E/3p(«)+&,/3V-^J'' ^i:/3p(»)+/3,/3*> -^^ 

3.2. Multivariate Interpolation preserving linear systems of PDEs 

Similar to the univariate case one can adapt the preceding algorithm to the inter- 
polation of multivariate functions, i.e. interpolate / by a polynomial p such that 
f =p, and 

(3.14) L,f{x) = L,p{x) for a; e 9. 

where Q = {xq, ■ • • ,xn} is the set of interpolation points, and the partial differ- 
ential operators are defined by 

(3.15) Uf{x) = J2 al,{x)d^f{x), = 1, • • • , fc. 

\a\<qi 

The procedure is analogue to that described in Section 2.2. (cf.also [7]). 

4. Approximation of global solutions of linear partial differential 
equations 

We refine the algorithm further in order to solve linear partial differential equa- 
tions globally. In this case the function u to be approximated is not known. In 
this section we shall simply describe an algorithm which constructs a polynomial 
which satifies a linear system of partial differential equations on an arbitrary set 
of interpolation points. It is not clear, however, if this polynomial approximation 
converges to the solution of the system. To ensure that and in order to estimate the 
rate of convergence we shall need the a priori estimates and regularity results. Note 
however, that the regularity constraints on the solution maybe low for problems on 
compact domains as any continuous solution functions u can be approximated by a 
families of polynomial functions approximating u. Therefore, principally, the fam- 
ilies of polynomial functions constructed here may approximate continuous global 
solutions in viscosity sense. An investigation of this problem will be considered 
elsewhere in a more general framework where we include some class of nonlinear 
problems. In order to make the basic ideas transparent we consider first scalar 
linear problems. We exemplify our algorithm first in the case of dimension n — 1 
and then generalize to the case n > 2. What we have in mind here are elliptic 
equations but we need the ellipticity condition only when we wan to prove that 
the family of polynomials construxted converges to the global solutions. Then we 
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exemplify our method in the case of a typical linear first order system. It is then 
clear how to generalize to systems of linear equations of any order. 

4.1. The case scalar second order equations of dimension n — 1 

We consider the simple boundary value problem 

d u du 

(4.1) Liu = a(x)-—r + b(x)-— + c{x)u = fix) on (d, e) C M, 

dx'^ dx 

with the boundary condition u{d) == Cd and u(e) — Ce (actually an ordinary differ- 
ential equation). If a{x) > A > for all a; S K, then we have an elliptic operator, 
but this is not an assumption which we need to construct an univariate polynomial 
which satisfies the boundary problem on the interpolation points. 

We start with the point d. We construct a hst of polynomial gmjTi > 0. We 
define the qm in substeps. Let po = oq. In order that po satisfies the boundary 
condition at a; = d we impose 

(4.2) Po = ao = Cd 
Next we define 

(4.3) pi{x) = ao + ai{x — d) 

In order to satisfy the second boundary condition we get 

c — c^ 

(4.4) pi{e) = ao + ai{e - d) = Cd + ai{e - d) = Ce ^ ai = — -. 

e — a 

It is clear that pi preserves the boundary conditions, i.e. p{d) = u{d) = Cd and 
p{e) = u{e) = Ce- Next let xq be the first interpolation point (any point in the 
interval (d, e). We want to ensure that 



(4.5) a(a;o)-7^(2;o) + Kxo)-jzi^^) + c{xo)p{xq) = /{xq). 



)— ^(xo) +b(xo)-—{ 
dx dx 

In order to ensure this, we define a polynomial which is an extension of po in three 
steps. First, define 

(4.6) P2{x) = flo + ai{x — d) + a4{x — xo)^(x — d){x — e) 
Plugging in and evaluating at x = xq we get 

(4.7) a{xo)2a4{xo - d){xo - e) + 6(a;o)ai + c(a;o)(ao + ai{xo - d)) = f{xo) 

Since oq, ai are known we get (recall that xq ^ d and xq ^ e) 

_ f{xo) - c{xo){ao + ai{xo - d)) - b{xo)ai 
2a{xo){xo - d){xo - e) 

Next define 

(4.9) p^ix) = P2{x) + a^ix - xo)(x - Xd)(x - x^). 
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Plugging in and evaluating at x ~ xq we get (assuming that ) 

■^iP3(a;o) = Lip2{xo) + a{xo)a3i2{xo - d) 
(4.10) 

+2(a;o - Xe)) + bixo)a3{xo - d){xo - e) = fixo). 

Hence, (provided that xq ^ d and xq 7^ e), 

(4.11) aa = .fM-LiP2M 

a(xo)(2(a;o - d) + 2(a;o - e)) + b{xo){xo - d){xo - e) 

Finally, finishing the first inductive step of recursive definition of the polynomial 
family {qm)meN 

(4.12) piix) ^ psix) +a2{x ~ d){x - e). 
Plugging in and evaluating at a; = xq we get (assuming that ) 

(4.13) LiP4{xo) = LiP3{xo) + 2a{xo)a2 + b{xo){{xo - d) + {xq - e)) = ,f(xo). 
Hence, (recall again that xq ^ d and xq ^ e), 

.--..N ^ fixg) - Lipsjxo) - b(xo){{xo - d) + {xq - e)) 

^ ' ' "' 2aixo){{x„-d) + 2{xo-e)) 

Now wc can define 

(4.15) qi{x)=p4{x) 

Next assume that the polynomials 91, • • • ,qk have been defined. This means that 
we have computed the polynomial coefficients gq, ai, • ■ ■ , a2+3fe. Then qk+i is de- 
fined via 

(4.16) qk+iix) = qkix) + (x - df{x - efn'[^f){x - xifzk{x), 

where Zk is a polynomial function which will be defined in three substeps. First, 
let 

(4.17) qk+i,i(x) = qk{x) + a2+^k+i){x - Xk+if{x - df{x - efU^^i-,{x - xif 

Plugging in leads to 

Liqk+i.ii.Xk+i) = Liqk{xk+i) + a{xk+i)2a2+3(k+i)ixk+i - d)^^ 
(4.18) 

{xk+i - efT^i=o{xk+i - xif = fixk+i). 

Hence, 

f{xk+i) - Liqk{xk+i) 



(4.19) a2+3(fe+i) = 



a{xk+i)2{xk+i - dY{xk+i - e)^n;*LQ(xfc+i - xi)^ 
Next, let 

(4.20) qk+i,2{x) = qk+i.i{x) + a2+3k+2{x - Xk+i){x - df{x - efn^^Q{x - xif 
We define 

(4.21) R{x) = (x - df(x - efUtoix - xif . 
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Plugging in leads to 

Liqk+i,2{xk+i) = -^i9fc+i,i(^fc+i)+ 
(4.22) 

a{xk+i)2a2+3k+2-£^Rixk+i) + b{xk+i)a2+3k+2-£^R{xk+i) = f{xk+i). 

Hence, 

/. ^g^ _ f{xk+i)- Liqk+i,i{xk+i) 

(i.Z6) a2+3k+2 - — . d2 „, r~"T7 X d nf T 

a[Xk+i)j7^R{xk+i) + b{xk+i)-^R{xk+i) 

Finally, let 

qk+i,3ix) = qk+i.2{x) + a2+3k+i{x - d)^{x - e)^n;*Lo(a; - xi)^ 
(4.24) 

= a2+3k+iRix) 

Plugging in leads to 

Liqk+i,3{xk+i) = Liqk+i,2{xk+i) + a{xk+i)a2+3k+i-^R{xk+i) 
(4.25) 

+b{xk+i)a2+3k+i-j^R{xk+i) + c{xk+i)a2+3k+iR{xk+i) = f{xk+i)- 
Hence, 

f. o«x f{xk+i) - Liqk+i^2{xk+i) 

(4.ZDJ a2+3k+2 — -^ -; ■ 

aixk+i)2^Rixk+i) + b{xk+i)-^R{xk+i) + c{xk+i)Rixk+i) 

It is clear how to proceed inductively in order to get a family of interpolation poly- 
nomials which satisfy the differential equation on an increasing set of interpolation 
points. Note, however, that we have not used any structural information about the 
coefficients at this point. This means that the equation may be ill-posed, and con- 
vergence cannot be guaranteed. 

4.2. The case of scalar linear partial differential equations 

For a positive integer k consider an equation of form 

(4.27) LkU = J2 ""(^)a^ = 5(a^)' 

\a\<k 



to be solved on the domain Q where 

(4.28) y 



on 



What we have in mind is an elliptic equation f order fc, but ellipticity is not re- 
quired in order to describe the algorithm which produces a family of multivariate 
polynomials which satisfy the equation on a set of interpolation points in fJ. El- 
lipticity becomes important when we want to show that the family of polynomial 
converges to the solution of the equation (assuming that there is an unique global 
solution). For simplicity of notation we consider the case fc = 2, i.e. the situation of 
p.9p . Assume that f ^ C^ and choose a discrete interpolation set 0f, C dfl. Then 
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we can apply the extended Newton algorithm of Seetion 3 in order to produce a 
polynomial pi, : M" -^ R such that 

Pb{x) — f{x) for all X e 06 
(4.29) 

|P| = |Pi for all a with lal < ? and x G Ob 

We assume that B^ — {xo6, • • • , XMb\ with Xi^ = {x\^^ • • • , x^^) and define 

(4.30) ^,{x)^TltLl,IiU^x^-xlf+\ 

Next let 6'i„f C fi \ Sfi be a set of interpolation points in the interior of fi. Let 

(4.31) Qint = {xo,- ■■ ,xn} ■ 

We enumerate (case k — 2) the q :— ^" „ -^'^ diffusion coefficients Oaj,--- , a^^ 
(arbitrary order), where we assume a/ — {oiii^ai2) and define first q polynomials 
Pq' ' (x), ^ = 1, • • • , (7. Let 

(4.32) pr'^i(x) =pb(x) + <I>b(a;)a„,(x"" -x^")(a;"- -x^-). 
Then we have 

(4.33) L2]5o*^(.To) = L2Pb{xo) + $h(a;o)(l + <5aiiai2)aai = /(a^o), 
which leads to 

_ f{xt)) - L2Pb(xo) 



(4.34) 



(f +<5aiiai2)*b(a;o) 



Having defined Pq' ' (x), • • • ,Pq ' (x) (and therefore computed a^j^, • • • , aa,) we 
define 

{4'35) p'^'»'+\x)=pt'''-\x)+M^W,,A^''^'^'^' -x'^,^'^^^^^ 

and evaluation leads to 

(. ofi^ „ /(xo)-£2Po'"''(p(a:o) 

Proceeding inductively we get app' '"^(x) which equals together with its derivatives 
up to order I the function / and such that the diffusion part of the operator applied 
to Pq' '"^(x) equals g at xq. It is now clear how this procedure can be extended such 
that an extended polynomial po (a;) equals together with its derivatives up to order 
I the function / and such that the total operator applied to ^0(2;) equals 5 at xq. 
As in Section 3 the ansatz for the interpolation polynomial pQ which satisfies the 
linear equation on the set of interpolation points O = {xo, • • • , x^r} then is 

N 

(4.37) pe{x) ^ Y. n}=in;'=i(^' ~ ^Ufp^(^)^ 

1=0 

where pi for i > 2 are then constructed as po above. 



18 Jorg Kampen 

4.3. The case of a linear hyperbolic equation 

We consider the hyperbolic equation mentioned above of the form 

(4.38) Lu = f in Q, 

where 

(4-39) Lu^Y. ^^^ 7^—1^. + E 7^ + ^(^)" 



] 



and (hij) is a symmetric matrix of signature (n, 1), if dimfi = n + 1. Note that the 
operator L can be transformed into the form 

(4.40) Lu=au + Liu, 

where Liu is some first order differential operator on Q. We assume the initial 
conditions 

(4.41) u ~ g and du = lo, 

where g and lo {1— form) are initial data. It is clear that the algorithm described 
in the preceding section can be used in the present situation. Later we shall see 
that energy estimates imply convergence of the scheme. 



5. Further refinements: collocation and parallelization 

Numerical experiments show that the coefficients of the recursively computed poly- 
nomials have to be computed with increasing accuracy in order to control effects 
of the truncation error of the coefficients of the polynomials. In the numerical ex- 
ample below, where we computed a polynomial approximation of degree 74 of the 
locally analytic function 

(5.1) ^-T^ 

1 + X 

and its derivatives up to order 3 on the interval [0, 5.4] such effects are not observed. 
However, if we increase the number of derivatives to be approximated up to order 
A: = 10 and increase the number of interpolation points, effects of truncation errors 
can be observed for polynomials of degrees larger than 200. The error increases 
as |a;| becomes large and truncation errors increase. This error can be reduced 
by a more precise representation of the computational approximation of the real 
numbers involved in the computation. However, as we point out in this section, 
we can compute m polynomials pj" ^ , • • • , p^™ of degree Ni,N2- ■ ■ Nm parallel 
which interpolate a given linear system of partial differential equations on some 
interpolation sets 0i,-- • ,&m using our basic algorithm, and then compute one 
polynomial py- q which interpolates the same linear system of partial differential 
equations on the set ^ O = 8i, U • • • , U8m- It turns out that this can be in such 
a way that the truncation error of the resulting polynomial py-- q is much smaller 
than in case of a direct extension of one polynomial pQ. using the basic algorithm. 
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We call this method the collocation extension of our basic algorithm. We shall 
assume that the sets of interpolation points are mutually disjunct, i.e. 

(5.2) e, n Oj = iff i 7^ j. 

It is clear that the computation of the polynomials p^ \ • • • , p^" can be done 
parallel and only the step of synthesizing has to be done non-parallel. Next we 
describe that step in case of two polynomials for simplicit of notation. Extension 
to TO > 2 polynomials will be clear from that description. So let Oi, 82 C fJ C M" 
be two discrete finite sets of interpolation points of a linear system of partial 
differential equations Lu = / to be solved on a domain fi and such that 0i n 82 = 
0. We write down the polynomial in the univariate case because this simplifies 
the notation, and the multivariate case is quite similar. Then we define a regular 
polynomial interpolation formula on 0i U Q2 by 

1^0=1 ^^Mj, (,«i_ «1).+i^^*=i (,«i_,«2).+iPei {X) 

\j k ' ^ J z -f 



(5.3) 






E, n,e{i,2} iM^ xf-f+w,^{x xfy 

--■■ Ej 962.61(2^)^61(2;) + hl{x) 



+ Ej ^ei ,62 i^)PB2 ix) + K {x) , 
where the constants a\j,,p G {1,2} are computed recursively as follows: For each 
j we can define Oq ^ = 0. If a{ p, • • • , al_^ are determined, then compute aj ^ via 

(5.4) E ( n Dlq^d,eS^,)Di-''PeAx^) ^ Dih^x,) 

l<r<l ^ ^ 

for each j. Note that this 'synthesis of polynomials' improves the computational 
power of our method dramatically. In the example below, where we approximate 
a simple locally analytic function 

(5.5) ^^T^ 

1 + X 
(with convergence radius 1) and its derivatives up to the third derivative on the 
interval [0, 5.4] with 19 interpolation points Oi = {fc0.3|fc = 0, • • • 18} we compute 
a polynomial of degree 74 in half a minute on a modest laptop machine. If we 
want to compute a polynomial which gives the same kind of approximation on the 
interval [0, 5836, 8] it will take several weeks. However, using parallelization and 
synthesis, and using the rough estimate that synthesis takes in average the same 
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time as building the 1024 basis polynomials of degree 74 on the intervals [0, 5.4] 
and [fc5.7, (fc + 1)5.7], A: = 1, • • • 1023 we need 10 steps of parallel synthesis of pairs 
of polynomials of cost of a less than a minute to get a regular approximation 
polynomial which is at least of degree 75776! It is clear fromthe preceding remarks 
how to extend this to the multivariate case (cf. also [7]). 



6. Convergence of polynomial approximations of global solutions of 
linear elliptic PDE and error estimates by a priori estimates 

Up to now we just considered (regular) polynomial interpolation on given sets of 
interpolation points. In this section we consider standard problems in the theory 
of linear partial differential equations and derive the convergence of our algorithm 
and error estimates (as the mesh size of the sets of interpolation points converges 
to zero). We start with elliptic equations and then consider hyperbolic problems. 
Similar results can be obtained for initial-value boundary problems for parabolic 
equations (since analogous error estimates can be obtained). In this case, however, 
it turns out that (at least for regular data) a WKB-expansion of the fundamental 
solution has better convergence properties and error estimates can be obtained by 
Safanov a priori estimates (cf. [8] and [6]). We shall consider application of our 
algorithm to this case in the next section. Note that Since to get an error from 
simple Taylor expansion in genera, because the interpolated function is unknown. 

6.1. Convergence for elliptic equations with regular data 

We consider the Dirichlct problem for elliptic equations, i.e. an equation of the 
form 

dx' 



(6.1) Lu=Y^ ««(^)7^ = /(2 

\a\<k 

on a domain Q, C R". coefficient functions 

(6.2) X -^ aa{x), 
and where u is given on the boundary, i.e. 

(6.3) ^\on=9- 

We consider the classical case where k — 2 and il is bounded. We assume 
uniform ellipticity, i.e. there exists a constant K > such that for all x G 17 

n 

(6.4) E«^^(^)^»^J-^^I^I'- 

In the classical case Schauder boundary estimates are available. We cite them 
in the context of a standard existence result, for a scalar function /i in fi we 
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introduce the norms 



k 



(6.5) l|/i|ir = EEll^''^llo 

J=0 \S\=j 

where 

(6.6) \\h\\o:=SMp\h{x)l 

and 

(6.7) \\h\\r+^ = \\h\\r + J2T.H'cf{D'h), 

where H^{J) is the Holder coefficient of a given fmiction / in il. We assume that 
the coefficient functions x — > aij{x) (diffusion terms), x — > biix) (drift terms), 
the potential term (x —> c{x)), and the right side x — > f{x) are uniformly Holder 
continuous (exponent a) such that 

(6.8) lla.,11^'^ < C, \\b.C < C, \\cC < C, WfC < C 
for some generic constant C. 

Theorem 6.1. Assume that conditions (j6.4[) and (j6.8p hold, and assume that c < 0. 
Furthermore, assume that dft belongs to C^^" and that g belongs to C\^d. Then 
the inequalities 

M\%^ <C(||g|| + h||„ + 11/11:;'^) 
(6.9) 

<C(|lg||+supa^|g|+Csupo I/I + 11/11^'^) 

hold. Furthermore there exist a unique solution u G C'2^^ to the Dirichlet problem. 

The interpolation polynomial pq described in the preceding section is by 
construction such that 

(6.10) L{u -Pe)^Yl «o(x) ^^"^~/^^ - A/(x), 

|a|<fc 

and 

(6.11) u-pe\g^^ Ag. 
It follows that 

Theorem 6.2. Assume the same conditions as in theorem 6.1.. Then 

(6.12) \\u - pe\\f+^ < C (II A,g|| + sup^o |A.g| + Csup^ |A/| + ||A/||^'^) 

Note that this implies an iy^-error even for the second derivatives of the global 
solution function, hence essentially an estimate in H'^{U). Even stronger results 
can be obtained if additional equations for the derivatives of u are considered (cf. 

H). 
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6.2. Convergence for a hyperbolic linear partial differential equations equation 

We consider again the hyperbolic equation mentioned above of the form 

(6.13) Lu ^ f on O C n, 

where 

ij ■^ i •' 

and (hij) is a symmetric matrix of signature (n, 1), if dimil = n + 1. We assume 
that some O C fi is bounded by two spacelike surfacesE^ and Eg and swept out 
by a family of spacelike surfaces Se(s). Recall that the initial conditions 

(6.15) u = g and du = lo. 

Let p be the interpolation polynom described above such that 

(6.16) L(m-p) = A/ onOc r2. 

(6.17) u — p ~ Ag and du — Auj. 
Then we use the following energy estimate 

Proposition 6.3. Let u solve the intial value problem (j6.13p . (|6.17p . Let 

(6.18) 0{s) ^O C] {t < s} 
(swept out by the spacelike surfaces llf.{s)). Then 

lois) \^\'dV < 
(6.19) 

/^.(^) \g\'dS + C{s - so) J^^ {\g\' + H') dS + C f^^^^ \f\'dV 

for s e [so,si]. 

This implies 

Theorem 6.4. With the same assumptions as in propostion 6.2. we have 

?{s) 



Iois)\^-P\'dV< 



(6.20) 

/^.(^) \Ag\^dS + C{s - so) J^^ (I A.g|2 + \cu\^) dS + C J^^^^ \Af\^dV 

for s e [so,si]. 

Hence the polynomial interpolation scheme described in Section 4 leads to 
L^-convergence. One can improve this scheme assuming regularity of solutions and 
considering systems of equations including equations for derivatives of the solution 
u (cf. H). 
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7. Applications to parabolic equations (connection to 
WKB-expansions) 

We summarize some results concerning WKB-expansions of parabolic equations 
(cf. [6] for details). Let us consider the parabolic diffusion operator 

where the diffusion coefficients aij and the first order coefficients bi in (|7.ip depend 
on the spatial variable x only. In the following let 6t — T — t, and let 

(7.2) {x,y) ^ d{x,y)>Q, {x,y) ~^ Ck{x,y), k>Q 

denote some smooth functions on the domain M" x R" . Then a set of (simplified) 
conditions sufficient for pointwise valid WKB-representations of the form 




(7.3) p{St, x,y)^ ,-— „ exp ^ + ^ c^ (a;, y)6t'' 



■§j^ — Lu = 0, with final value 



for the solution (t, x) — > p{St, x, y) 

(7.4) 

u[0,x,y) = 5{x-y), 

is given by 

(A) The operator L is uniformly elliptic in R", i.e. the matrix norm of (a^ (x)) 
is bounded below and above byO<A<A<oo uniformly in x, 

(B) the smooth functions x — > aij{x) and x — > bi{x) and all their derivatives 
are bounded. 

For more subtle (and partially weaker conditions) we refer to [6] . We consider the 
case where there exists a global transformation to the Laplace operator. If we add 
the uniform boundedness condition 

(C) there exists a constant c such that for each multiindex a and for all I < 

i,j,k < n, 



(7.5) -77^, t; — <cexp(c|x 
^ ^ ax" dx°' ~ ^ 

then the function cP = (x — y)^ (in the transformed coordinates and Ck equals its 
Taylor expansion around y E R", i.e Ck,k > have the power series representations 

(7.6) Ck{x,y) =J2a'^k,a{y)Sx°',k>Q. 
Moreover Ck,k > are determined by the recursive equations 



(7.7) -!l + Ud' + lj2ij2{a..,ix) + a 



^X, \ dcQ 
dxi 



3^i^))^ I 7rr(2^'y) 
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where the boundary condition 

1 



(7.8) CO (y, y) - - - In y/det{a%j)) 

determines cq uniquely for each y G M", and for fc + 1 > 1 we have 

(fc + lK.+i(x,y) + iE.,«..(-)(4-%f + %%t') 
(7.9) 

with boundary conditions 

(7.10) Ck+i{x,y)^ Rk{y,y) a X ^y, 

Rk being the right side of (|7.9p . In case Uij — Sij we have the representations 

(7.11) d^(a:,^/)=^(x.-2/.)^ 



(7.12) co{x,y) ^y^iyt-x,) 6, 

"'0 



{y + s{x-y))ds, 



and 

(7.13) Ck+i{x,y)= Rk{y + s{x-y),y)s''ds, 

Jo 

Rk being again the right-hand-side of (|7.9p . The integrals can be taken out if the 

functions x —^ bi(x) are given by multivariate power series and error estimates 

for the truncation error in space and time are obtained (cf. [SI [5]). However, even 

if the coefficient functions are analytic, i.e. equal locally a power series, it is not 

possible to approximate such a function globally by their Taylor polynomial. As 

an example consider the equation 

/r, , .X du 1 . v^ 1 du ^ 

(7.14 A-u-> T^=0 

^ ' dt 2 '^ l + xidxi 

Here, the coefficient functions 

(7.15) x,~>-^^b,{x) 

i \ Xi 

are univariate locally analytic function with convergence radius 1. Such type of 
equations occur in praxis of finance (cf. [IJIS]). In order to obtain an approximation 
of the WKB-expansion say up to order 5, i.e. compute the coefhcient functions 

(7.16) cc ^ Cfc(x,2/),fc = 0, • •• ,5, 

we need a global approximation of the functions (|7.15p and their derivatives up 
to order 10! This is due to the recursion equations for the Ck,k > 1 which involve 
second derivatives of Ck-i- If we have 20 interpolation points on the 2:-axis this 
implies that our regular interpolation algorithm computes a polynomial of order 
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231. We do the computation in a more modest example in order to keep the 
resulting polynomial representable on one page in the following section. 



8. A numerical example 

The following polynomial is a similtaneous approximation of the function 



/: [0,5,4] C 
(8.1) 



and its first, second, and third derivative on the domain [0,5,4] with 19 interpo- 
lation points. Hence the degree of this univariate polynomial is 74. Note that the 
convergence radius of / is 1 . 



75 

(8.2) p^^{x) = Y. «"(^ - :cmdiv4)"""°''^nj2.t;^4-i(^ _ ^^Y 

m=0 



Note that 



(8.3) 



dx" V 1 + a; 



1 ^ I (-l)"n! , , 1^„ , 



This leads to the values ag = 1, ai = — 1, a2 = 1, and 03 = —1 for the coefficients 
of our interpolation polynomial at xq — 0. Note that the coefficients a^ of the 
interpolation polynomial tend to become smaller for large indexes i as you would 
expect. 
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ao = 1.0 

ai = -1.000000000000, aa = 1.000000000000, as = -1.000000000000 

a4 = 0.769230765432, a^ ^ -0.591715921811, ag == 0.455165657066 



0-7 

aw 
ai3 
ai6 



= -0.350124490177, ag = 0.218822618520, ag = -0.136753442090 
= 0.085452696109, an = -0.053404312398, ai2 = 0.028144617397 
= -0.014953338935, au = 0.008262188243, ais = -0.005370784216 
= 0.003734873988, air = -0.003633027430, aig = 0.004645502211 
ai9 = -0.006813570816, a2o = 0.007086610952, a2i = -0.007144432312 
a22 = 0.006238342564, aas = -0.002646146059, 024 = -0.002374282360 
025 = 0.008387675067, a26 = -0.015766978592, 027 = 0.024857498610 
== -0.025373687351, 029 = 0.025025735340, ago = -0.023974098174 
1 ^ 0.022321168853, a32 = -0.015945627926, 033 = 0.011155207224 
= -0.007619506803, 035 = 0.005069120726, 035 = -0.002759684498 
.001479716734, 038 = -0.000790172686, 039 = 0.000430223475 
0.000208511304, 041 = 0.000106279314, 042 = -0.000056281013 
.000028862889, a44 = -0.000011153733, 045 = 0.000002201139 
a46 = 0.000002233629, a47 = -0.000004202246, a^s = 0.000003699371 
= -0.000002870941, a^o = 0.000002068390, 051 = -0.000001402599 
= 0.000000753699, a53 = -0.000000375935, 054 = 0.000000159621 
055 = -0.000000037499, age = -0.000000015690, agr = 0.000000032004 
= -0.000000031616, a^g = 0.000000023406, ago = -0.000000010968 
1 = 0.000000001564, ae2 = 0.000000005590, a63 = -0.000000011521 
a64 = 0.000000012190, ags = -0.000000012095, age = 0.000000011769 
aes7 = -0.000000011467, aes = 0.000000008698, aeg = -0.000000006652 
== 0.000000005132, an = -0.000000003988, 072 = 0.000000002523 
= -0.000000001600, a74 = 0.000000001015, 075 = -0.000000000643 



028 

03 

034 

(8.4) 

a37 = 

040 = - 

043 = 



049 
a52 



058 

ag 



070 
073 
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9. Conclusion 



We have designed regular polynomial interpolation algorithms and variations which 
produce families of multivariate polynomials which solve linear systems of partial 
differential equations on arbitrary sets of interpolation points. In our basic al- 
gorithm the members of the family of polynomials are defined recursively each 
being an extension of the preceding member in the sense that the preceding mem- 
ber agrees with a given member on the set of interpolation points on which the 
preceding member satisfies the linear system of partial differential equations. We 
have shown that the family of multivariate polynomials has the global solution as 
its natural limit if some a priori information on the system of partial differential 
equations is available. The information needed can variate from case to case. In 
any case a solution should exist. We have shown how to use a priori estimates 
of elliptic equations and of hyperbolic systems of equations in order to obtain 
error estimates adapted to the regularity of the solution. Similar is true for par- 
abolic equations. All this makes our approach compatible with new techniques 
like sparse grids or weighted Monte-Carlo algorithms developed in order to treat 
systems of higher dimension. In case of parabolic equations we showed how regu- 
lar polynomial interpolation of known functions can be used in order to compute 
higher order approximations of WKB-expansions of fundamental solutions. We 
also constructed extensions where the algorithm is parallelized on different set of 
interpolation points an showed how these partial polynomial approximations can 
be patched together to one multivariate polynom which fits the given system of 
linear partial differential equations on the union of sets of interpolation points. 
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